use "${clean}analysis.dta", clear
keep if north == 1 
drop if congress > 55
drop if state == "District of Columbia"

//standardizing variables
egen birth_state_surban = std(birth_state_urban_percent)
egen birth_state_srugged = std(birth_state_ruggedness)
egen birth_state_spostoffice = std(birth_state_postoffice_area)
egen sborn = std(born)

// creating a variable to keep the sample constant
cap drop sample
gen sample = 1
foreach var of varlist $state_bcontrols $life_controls {
	replace sample = . if `var' == .
}
keep if sample == 1

// dropping third parties
drop if partyid == 2 | partyid == .

// dropping WHIGS
replace partyid = . if party_name == "Whig"
replace partyid = 1 if bioguide == "B000376"
replace partyid = 1 if bioguide == "S000261"
replace partyid = 1 if bioguide == "S000801"


// creating some variables
local c = 1860

gen period = year > `c'
replace year = year - `c'


// matrices to store the results for graphing
matrix define I = J(2, 3, 0)
matrix define S = J(2, 3, 0)



* regressions
* ------------

// republican baseline mixed model
mixed singular c.year##i.period##ib1.partyid ${state_bcontrols} ${life_controls} if sample == 1 || _all:R.congress || bioguide_id:

qui lincom 1.period // intercept change at cutpoint
matrix I[1,1] = `r(estimate)'
matrix I[1,2] = `r(estimate)' + 1.96*`r(se)'
matrix I[1,3] = `r(estimate)' - 1.96*`r(se)'

local beta_r = `r(estimate)'
local beta_r : display %4.3f `beta_r'
local se_r = `r(se)'
local se_r : display %4.3f `se_r'


qui lincom 1.period#c.year // slope change at cutpoint
matrix S[1,1] = `r(estimate)'
matrix S[1,2] = `r(estimate)' + 1.96*`r(se)'
matrix S[1,3] = `r(estimate)' - 1.96*`r(se)'
	




// democratic baseline mixed model
mixed singular c.year##i.period##i.partyid ${state_bcontrols} ${life_controls} if sample == 1 || _all:R.congress || bioguide_id:

qui lincom 1.period // intercept change at cutpoint
matrix I[2,1] = `r(estimate)'
matrix I[2,2] = `r(estimate)' + 1.96*`r(se)'
matrix I[2,3] = `r(estimate)' - 1.96*`r(se)'
	
qui lincom 1.period#c.year // slope change at cutpoint
matrix S[2,1] = `r(estimate)'
matrix S[2,2] = `r(estimate)' + 1.96*`r(se)'
matrix S[2,3] = `r(estimate)' - 1.96*`r(se)'









* graphing 
* --------
svmat I
svmat S


// temp variable to index estimates
cap drop temp*
gen tempn = _n if S1 != .

gen temprep = tempn == 1


twoway ///
	(scatter I1 tempn if temprep == 1,	msize(large) mlwidth(none) mcolor("${orange}") 	msym(D) yaxis(1 2) ) ///
	(scatter I1 tempn if temprep == 0,	msize(large) mlwidth(none) mcolor("${green}") 	msym(D) yaxis(1 2) ) ///
 	(rcap I2 I3 tempn if temprep == 1,   lwidth(thick) lcolor("${orange}") 	lpattern(solid) yaxis(1 2)) ///
 	(rcap I2 I3 tempn if temprep == 0,   lwidth(thick) lcolor("${green}") 	lpattern(solid) yaxis(1 2)) ///
	, ///
	text(-.05 1 "{&beta} = `beta_r'") ///
	text(-.125 1 "se = `se_r'") ///	
	yline(0, lcolor(gs2) lpattern(dash) lwidth(thick)) ///
	legend(off) ///
	xlab(1 "Republicans" 2 "Democrats") ///
	xtitle("") 	///
	ylab(none, axis(2)) ///
	ylab(-.45(.15).45, axis(1) grid glcolor(gs7%50)) ///
	ytitle("") ///
	xscale(range(0.5 2.5)) ///
	subtitle("Intercept Shift: 1860 - 1861", box bexpand bcolor(gs2) color(white) size(medsmall)) ///
	graphregion(margin(vsmall)) ///
	name(intercepts, replace) nodraw
	
twoway ///
	(scatter S1 tempn if temprep == 1,	msize(large) mlwidth(none) mcolor("${orange}") 	msym(D) yaxis(1 2) ) ///
	(scatter S1 tempn if temprep == 0,	msize(large) mlwidth(none) mcolor("${green}") 	msym(D) yaxis(1 2) ) ///
 	(rcap S2 S3 tempn if temprep == 1,   lwidth(thick) lcolor("${orange}") 	lpattern(solid) yaxis(1 2)) ///
 	(rcap S2 S3 tempn if temprep == 0,   lwidth(thick) lcolor("${green}") 	lpattern(solid) yaxis(1 2)) ///
	, ///
	yline(0, lcolor(gs2) lpattern(dash) lwidth(thick)) ///
	legend(off) ///
	xlab(1 "Republicans" 2 "Democrats") ///
	xtitle("") 	///
	ylab(none, axis(1)) ///
	ylab(-.06(.02).06, axis(2) grid glcolor(gs7%50)) ///
	ytitle("") ///
	xscale(range(0.5 2.5)) ///
	subtitle("Slope Change after 1860", box bexpand bcolor(gs2) color(white) size(medsmall)) ///
	graphregion(margin(vsmall)) ///
	name(slopes, replace) nodraw
	
	
graph combine intercepts slopes, scale(1.5) xsize(100) ysize(50) graphregion(margin(r=-10 l=-10))
graph export "${output}FigD9_congress_north_demrep_statecontrols.pdf", as(pdf) replace

	

	
	
	







